Dynamics of inertial particles in a random flow with strong permanent shear 
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We consider advection of small inertial particles by a random fluid flow with a strong steady shear 
component. It is known that inertial particles suspended in a random flow can exhibit clusterization 
even if the flow is incompressible. We study this phenomenon through statistical characteristics of 
a separation vector between two particles. As usual in a random flow, moments of distance between 
particles grow exponentially. We calculate the rates of this growth using the saddle-point approxi- 
mation in the path-integral formalism. We also calculate correction to the Lyapunov exponent due 
to small inertia by a perturbation theory expansion. 
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I. INTRODUCTION 



Let us consider a system of small identical inertial particles suspended in a random fluid flow. The flow consists 
of random fluctuations of a fluid velocity on the background of a strong permanent shear. This system can serve 
I""!' as a model of various natural and experimental situations: plankton colonies in the ocean [l| or water droplets in 
\ warm clouds (which is useful for estimating rain initiation times 0,111). The most dramatic common feature of those 
1 situations is clusterization of inertial particles advected by a random flow. 

At first sight, one would expect an incompressible random flow to mix advected particles. However, due to inertia 
particles' velocity field differs from the fluid velocity field and becomes compressible. This compressibility leads to 
" unmixing" : emergence of clusters and voids in an initially uniform distribution of particles (T], 0] . 

Up to now, effect of inertia on a dynamics of advected particles was studied for isotropic fluctuating flows. For 
example, aggregation of inertial particles in a random compressible flow in one, two and three dimensions was consid- 
ered in '5'], [g'] and [tI respectively. It was shown that small inertia reduces mixing, that means it produces negative 
corrections to The Lyapunov exponents, whereas large inertia can turn the largest Lyapunov exponent to zero. This 
I is interpreted as a phase transition from aggregating to non-aggregating phase. It should be mentioned that since 
\^ . the velocity field of particles deviates from the fluid velocity field in the presence of inertia, a sum of the Lyapunov 
04 ■ exponents is generally non-zero. Moreover, this sum is always non-positive and can be interpreted as minus the rate 
of the mean entropy production per unit volume Q . The statistics of the concentration fluctuations was studied 
■ in |10| . In (TT| the corresponding multifractal set in phase space was studied using concepts from dissipative dynam- 
ical systems theory (as opposed to behavoir of passive tracers, which can be thought of as a conservative dynamical 
system). 

However, in many experimental and natural systems fluctuating flow does not exist alone, but against the back- 
ground of a strong shear that changes slowly. Processes we are interested in take place at scales much smaller than 
the viscous scale of turbulence, so the velocity field is spatially smooth. We model it by a permanent shear fiow with 
superposed weak gaussian spatially-uniform fluctuations of velocity gradients. Fluctuating part of the flow velocity is 
thereby short-correlated in time and linear in space. 

Dynamics of tracers (particles without inertia) in a shear flow combined with weak fluctuations was studied in 
hj,hj. The pure shear flow doesn't produce exponential growth of the distance between tracers, the distance grows 
only linearly. However, when combined with weak fluctuations of a flow velocity, the separation between the particles 
exhibits a non-trivial dynamics. 

Consider two particles suspended in such a flow. Because of the shear component of the flow one of them is moving 
on average faster than another (in this paragraph we project all the motion on the direction of the shear velocity) . As 
the shear flow makes the first particle move faster, most of the time it is situated to the right of the second particle. 
However, from time to time noise can overtake the average velocity difference and place the second particle to the 
right of the first. This situation will persist for the short period of time, after which the first particle will again outrun 
the second. For separation vector this means the following behavior: most of the time it spends in the region of small 
angles of order (j)Q — {D/ sY^^ along the shear direction, where s is shear rate and D is intensity of fiuctuations. From 
time to time, with an average period T ~ (_Ds^)~^/'^ it tumbles, making a half-turn and again for a long time it stays 
directed along the shear velocity. 

This dynamics produces a stationary angular PDF that is mainly localized at small positive angles. However, due 
to tumbling, it also has an algebraic tail 7^((/>) oc (f>~'^ dX <j)Q <^ (j) <^ 1. The Lyapunov exponent is proportional to 
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{Ds'^y^^ and n-th moment of distance grows exponentially with the rate {n'^Ds'^y^^ for n ^ 1 p^ . 

The task of this paper is to study a system, in which both effects are present: particles have small inertia, and 
the two-dimensional random flow has a strong permanent shear component. In our model there are four parameters 
that have dimension of time: the inverse fluctuations intensity D~^, which is much large then the inverse shear rate 
s^^, which is in turn much larger than the measure of particle inertia r. The smallest is the correlation time of 
fluctuations, which we completely ignore in this paper, treating fluctuations as a white noise. 

In Section [II] we calculate a correction to the Lyapunov exponent due to particles' inertia. We derive the system 
of stochastic equations, from which the Lyapunov exponent is expressed as an average of one of the variables. The 
corresponding Fokker-Planck equation is then transformed to a Schrodinger equation of a harmonic oscillator with 
a non-hermitian perturbation. We then perform perturbation theory calculations. The correction to the Lyapunov 
exponent turns out to be proportional to — tA^, where A ^ (Ds^Y^^ is the Lyapunov exponent of the tracer. In Section 
lllll we calculate the growth rates A„ of high moments of a separation between two particles. For large moment numbers 
this can be done analytically in a saddle-point approximation using path-integral formalism |l5| . We reproduce the 
result from inertialess case A„ = {n^Ds'^y^^ of [l3| and find the correction to it due to a small inertia, which is 
~2nT{nDs^)^/^ . Finally, in Section ITVl we summarize the results and discuss possible paths of a further work. 



II. LYAPUNOV EXPONENT 



The Lyapunov exponent is defined as the limit A — lim t ^ \nR{t). It can be expressed from the rates of exponential 



growth of the distance, A„ = lim t In (i?"), at small n: 
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A-^|„=o. (2.1) 
an 

For tracers this quantity was calculated in [isj . In this section we calculate the first correction to it due to a small 
inertia. 

Let us describe here a formal model, in the framework of which we'll work further. Velocity field, in which particles 
are suspended, consists of two components: permanent shear flow along a:-axis 

Vx = sy (2.2) 

and fluctuating part Vi — aijRj. So the full tensor of the velocity gradients is as follows: 

Y.ji{t) = sSjxSiy + (Tji{t). (2.3) 

It should satisfy the incompressibility condition TrE = 0. 

Fluctuating part of the velocity is assumed delta-correlated in time, isotropic and incompressible 

{(Tik{tl)(Jjn{t2))^D{?i6ij6kn-hk6jn~~hn6]k)6{ti-t2), (2.4) 

It is worthy of notice that in a shear flow the only relevant component of a is a2i, so the isotropic tensor structure of 
(|2.4p is chosen only for convenience The fluctuations are assumed weak relative to the shear, that is D ^ s. 

Spherical particles in an incompressible flow u(r,<) are subject to Stock's drag force f = 67r?7a(u(r) — r). Consider 
two such particles with a separation R(t) between them. Velocity field is smooth in our model, so in the first 
approximation u(r -(- R) — u(r) — SR. This leads to the following motion equation for R 

rii + R^SR. (2.5) 

Here r = m / [Q'K'qa) is a measure of particle's inertia, which we assume small, namely t <^ D. It has a dimension 
of time and can be interpreted as the "response time" of a particle. We assume that all particles are identical and 
don't interact with each other. 

We pass to new variables: an angle </> respective to the shear direction, p — \nR and Yi — p^Yi — 4>. Equation (|2.5|) 
transforms to 



tYi = -Yi + t{Y^ - Y]2) + s sin cos </> -I- 
TI2 = -Y2 - 2tYiY2 - s sin^ + ^2 

= ^2 



(2.6) 
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Though in (|2.5|) noise was multipUcative, in (|2.6I) it has become additive: ^1,^2 are independent white noises that are 
hnear combinations of elements of a (the price we pay for it is that the equations become nonhnear). As long as the 
dynamics of 0, Yi, 1^2 is governed by the system of stochastic equations (|2.6|) . their joint PDF satisfies Fokker-Planck 
equation r9tP((/), Yi, 1^27 ~ LP{(f>, Yi,Y2, t). The Fokker-Plank operator L is 

L = ~TY2d^ + di (yi - s sin cos + ^di^ + 82 (y2 + s sin^ + + tOi {Y^ - K^') + 2x^92^2 (2.7) 

P(0, Yi,l2) finally takes a stationary form, while p grows infinitely with a stationary stochastic increment. Thus we 
look for zero modes of L. 

In order to make perturbation theory calculations more convenient, we pass from the Fokkcr-Planck to a Shrodinger 
equation using a substitution 

P(0,xi,X2) = Q(0,a;i,a;2)exp(-a;?/4-x^/4) , (2.8) 



where xi = \Jt jD {Y\ — s sin cos (jj) , X2 = \Jt/D {Y2 + s sin^ 0) . 
Operator L transforms as follows: 

H = _e(^i+^')/4Le~(^i+^^)/4 = _g [43,^ _ ^1/2 - xixl/2 - X2d^ + 22;ia;2a2 + {xl - xl)di^ 



-a 



((4 -xl + 2xidiY^ - ^ + 02 _ 202a;ia2 + X2di + 202^291 



A 



2_|_ 2 



where e = V-Dr, a — s^Jt/D. We assume that e <C 1, ea ^ 1. 

The resulting hamiltonian H corresponds to a two-dimensinal harmonic oscillator with a non-hermitian perturba- 
tion. We express it through oscillator creation-annihilation operators a ~ x/2 + d^, a) ~ x/2 ~ dx 



'H(°'>=alai 



+ 0302 

-X2dcj, - 



4 a? 



Hi^) ^ -ea 2 sin 20 02 



2a[ ai — a[ 

t ■ 9 
2aia2 sin ( 



I ,-,;„2 



i| ) — 2a2ai(a2 -I- 03) oc -y/r 
— a|a2 cos 20 — xioj sin 20 oc t 



(2.9) 



Below in this section we will find corrections to a ground state wave function by perturbation theory expansion in 
small parameter e. Let us notice that our perturbation theory is degenerate, because acts as identity on the 
space of functions of and thus its eigenfunctions can have an arbitrary angular dependance. 
We look for a new ground state as a sum over non-perturbed oscillator eigenfunctions 



Q{X1,X2,<P) = ^/n(0)|? 



(2.10) 



Further we denote by fi%2 a correction of k*^ order to /(0) corresponding to the oscillator eigenstate with quantum 
numbers ni,n2. Obviously foo ~ 1, and all other fi are small at least as 0(e). As we will see below, an equation 
that determines the proper zero-approximation wave function coincides with stationary Fokker-Planck equation for 
inertialess problem. 

We substitute series (|2.10p into HQ = and obtain the following system of equations 



0, 



(0) 



ki 



(2.11) 



Here we denote Hnk = {nin2\H\kik2) - matrix elements of H between the non-perturbed oscillator states (they 
still act as operators on the space of functions of variable 0). 

We expand (|2.1ip successively in e and obtain the following corrections to /„ for ni + n2 > 



' f(l) 

Jn — 


^^io' f(0) 
_E„ J 00 ' 






/ rr(2) 


^ ^nk "kO 














nEk 



AO) 
J 00 



,(1)h-(1) 




(2.12) 



H. 



(2) 
00 
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Here and in what follows we sums n, fc, m are over all oscillator states except the ground state. Let us notice that 
^(2) preserves parity of sum of oscillator quantum numbers 7ii + n2, and i?^^-* inverts it, so /„ with even ni + 71,2 have 
corrections only of even order in ^/r, and vice versa. 

It's a bit more difficult to calculate /oo- Expending (|2.1ip up to the 2nd order, we get an equation for /qq'': 

rr(2) AO) _ -"ofc ^ kO AO) /r, , o\ 

-"00 JOQ ~~ j^^ JOO ■ l^-^-Jj 

Non-zero matrix elements entering this expression are (00|i?|00) ~ —ad^sva^cj) and (00|-ff|OI) = (01|-ff|00) = —ed^. 
Thus (j2.13p coincides with inertialess Fokker-Plank equation for an angular PDF 

a^(ssin2 + i?a^)/^°^ =0 (2.14) 

This equation was analyzed in [l3||. In the region of small angles its normalized solution is {s/DY^'^F ((s/Z?)^/'^0), 
where 

00 

F{^) = cj exp (-T^VS - (C - ^)V3) , C = ^-^^^^ (2.15) 


With this PDF, one can show that the inertialess Lyapunov exponent A = s((/)) is 

(2.16) 

r(i/6) ^ ^ 

In order to calculate a correction to this exponent, we have to expand (|2.1ip to the 4th order. Corrections of odd 
order to /oo vanish, the first correction to /oo is of the second order and satisfies the equation 

/ fj'-^^ fj'^^'^ fj'-^^ J- fj'-^^ H^'^'^ fj'-'^^ ^ fj'-'^^ H^'^'^ fj'-^^ 77-(l) 77-(l) H-(l) 

tt(2) r{2) _ _ I \ " J^On J^nm-nmg + J^Qn J^nm-nmq + J^Qn -"rim-"mo \ ' -"On ^nniJ^^k^kO 1 f{0) icy -, 

^00 JOO — 2^ p P p p F MOO l^-J^'J 

\n,m n,m,k I 

Operator in the right-hand side is 

- [Ds {-84, sm2(j>d4, + 84, cos 20 - 5^ sin 09^ sin^S^) - (29^ + 9^)) (2.18) 

The term with second derivative renormalizes the fluctuations power D : D' = D{1 + 2Dt). 

We expand p.l7p in small (f> and search for a solution in the form /qq"* ~ {s/ DY/^F^'^^ ((s/I?)^/'^0) . In the main 
order in D / s: 

{drjv^ + d^) f(^\ti) = -T{Ds^f'^ {-dr, + d,,7jdr,vd^ + 2d„vd„ + 8*) F^^Ht^) (2.19) 
We look for a solution of (|2.19p that doesn't violate normalization condition, that means / d(f>F<-^\<j)) = 0. Its 



explicit form is 

f'^^\t]) = T(Z?s2)i/3g-.rV3 J J dT^^e-'^l/^ J ^7736"^'/^ (Mt^i) ~ /i(773)) , (2.20) 

— OO — oo — oo 

where h{r]) = (l - Sr/dn - rj^d^ - d^) F^^^tj) 

Our aim is the Lyapunov exponent, which is given by 



A = (n) = VD/^i^i) + s('^> (2-21) 

One can see that the correction to The Lyapunov exponent comes from two sources: correction to the mean angle 
and non-zero (xi). The correction to the mean angle can be calculated from the correction to the angular PDF 

2 2 '^^ 
(0) = y dxidx2d(j>cj)exp (^-^1±^^ Q(si,X2,<^) = I 0/oo(0)d'/' = (0)o + J c^f^l\cl^)d<P (2.22) 



5 



+ 00 f] +00 

Ci = - / dT^T^e-"'/^ [ driie'^'/^ [ drne'"'/^ I drj^e"'^/^ {h{rji) - h{r]3)) « 9.8 (2.23) 



So, increment correction due to change of the mean angle is SXi — s6<j) — — Cir(Z?s^)^/'^ 

The second increment correction is related to the fact that averaged Yi is no longer equal to s((^) since we take 
inertia into account 



22 ^"^ 
{xi) = J dxidx2d<j)Xiexp ^- ^^ 4^ Q{xi,X2,(f>) = J 





Qixi,X2,(t>) = / .fio{(t>)d(t> (2.24) 



Full derivatives vanish after an integration, so the first order correction /Iq"* ~ d^f^ doesn't contribute to the answer. 
The second-order correction to /lo is zero, consider the third order (I2.12p : 



fw = - ( e'aa^ 1 + ^ cos 20 j + e^a sin 20 + ^e^S^ ) /^^ (2.25) 



2ir 

{xi) = -e^a j d0sin20/^°V) (2-26) 


We see that correction to the Lyapunov exponent due to non-zero (a;i) is much smaller than the correction due to 
change of the mean angle and can be neglected as long as D/s 1: 

5A2 = \ —{xi) -2DsT{(j)) = -2DT{Ds^f'^ < 5Xi. (2.27) 

V T 



Thus we have calculated the correction to the first Lyapunov exponent. The correction is negative, so inertia reduces 

WW 



mixing like in an isotropic situation. Using result for inertialess Lyapunov exponent A = ^i/L D^/^s^/^ from [T 



we can represent the answer in the relative form 

5X/X = - ^'^'i^^^ Ar -46.5Ar (2.28) 

III. LARGE N SADDLE-POINT APPROXIMATION 

In the previous section we have calculated the Lyapunov exponent, which is determined by the growth rates A„ of 
{R"{t)) at small n (j2.ip . and in this section we'll calculate growth rates for large n. The method we exploit to do 
this is based on the Optimal Fluctuation principle (saddle-point approximation in a functional space) . Its application 
to problems of statistical hydrodynamics was developed in papers [Tsl - fTTj . To make the saddle-point approximation 
more transparent, we use path-integral formalism for averaging over trajectories of stochastic processes. 

First we develop a path-integral representation of i?^" averaged over the velocity statistics. Random velocity enters 
equations (12. 6p through two independent white-noise stochastic processes, so we can express the average as 

(i?2"(t)) - J DCMR^'^VlCiMCi] cx e^-* (3.1) 

where P is a probability measure of a white noise 

m = ^^v(-^j dt'eY (3.2) 
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We express the noise from Langevin equations (|2.6p and the functional integral over the noise transforms to the 
functional integral over 0, Yi 



J DY,D<j>exp{-S[Y,,cl>]) 
with an action determined by the original Langevin equations 



(3.3) 



S[Yi,cf,] = J dt' 





tYi -Yi+T{(j)^-Y^^) + s sin (j) cos ( 



2TYi(j) - ssin 



2nYi 



(3.4) 



At large times and n 3> 1 the main contribution to this integral is given by the saddle-point approximation 15]. This 
means that we calculate optimal trajectory (j)*(t),Yi(t),Y2{t) which minimizes 5, then in the main approximation 
(I3.3P equals to exp {—S[(j)*, Y^, 1^2*])- ^.t first we will follow this procedure for slightly more general form of Langevin 
equations than (|2.6p in order to use it for both inertial and inertialess cases. 

Consider the following system of stochastic equations 



|0 = /2(0)+6, 

where ^1,^2 are independent white noises. The corresponding action is 



(3.5) 



t t 



- 2nYi 



(3.6) 



Motion equations and boundary conditions are obtained by varying this action with respect to 5Yi , 

■Yi-h = 2nD, 

^(0-/2(0)) |o,t = 0. 
The second equation describes Newtonian motion in the potential 

f/(</)) = 2nDh{4>) - fl{c^)/2, 

thus we can introduce conserved "energy" 

[^f/2 + U[<P)^E. 
Now we have solution for optimal trajectory in an implicit form t = J 



(3.7) 

(3.8) 
(3.9) 

and substitute it into the action 



5*: 



I dt' (^-2l5 (/2(0) - V2E-4nDMcb) + fi{cf,)y + 2Dn^ + 2n/ij - 

= ildt' (/!(</>) - 4n^/i (</>)) - 2DnH - i/2(0)0|^ + Et/D 


Using energy conservation relation, the first term can be rewritten as 2U = 2E — (jp. Finally the growth rates are 
expressed as follows 



X2nt f + -^/2(</')0lo + Et/D + 2DnH 



(3.10) 



Below we apply this formula to concrete /i, /2 corresponding to inertialess and inertial dynamics. 
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In the inertialess case (r = 0) from (|2.6p one can see that fi{4>) = ssin0cos0, /2((/)) = —ssiin^ (j). The potential 
U{(/)) — nDssin2(/) — sin** 0/2 has a maximum at = (nD / s)^/'^ . Boundary conditions require that 

TP 

sin 20(0) = sin 20m = — — , (3.11) 
nUs 

so the optimal trajectory is tumbling from some initial positive 0(0) to some angle 0(t) = kn + 0(0). As long as we 
keep time t is fixed, action doesn't depend on the number of turns k. We investigate behavior of the system at large 
times, so optimal trajectories which we consider should take long time. This is achieved if a gap between the energy 
and the maximum of the potential is small : t ^ \n{E — U{(j)m)), so E ^ Umax = ^{n^D^s'^y^^ . 

The first term in the action p.lOp doesn't grow with time, because it can be estimated from above: 
/ d(j)y^2{E — U{(j>)) < nk ■ 2maxAJ7(0). So we neglect it as compared with Xt. The second term in p.lOp van- 
ishes due to boundary conditions p. lip . Finally, increment is determined by the energy alone 

X2n^ E/D=^{n'^Ds^)^^^ + 2n^D. (3.12) 

Consider now a correction linear in r. Formally, order of the motion equation increases, but all terms with high 
derivatives are multiplied by small parameter, so we can substitute all higher derivatives from the inertialess solution 



- ^2{E-Uom, 4> - -C/o(0) (3.13) 

and prevent the order of the motion equation from increasing. Increment is still determined by an energy, which is 
very close to maximum of the potential. 

C/(0™) = C/o(0l'!> + <50) + ,Si7(0(^) + 50) = f/o(0„) + <5f/(0W) + 0(t2) (3.14) 

From dlH) and (gH) follows 

6U ^2nDS,h- f26f2, 

(5/2 = -r(2yi0 + 0). (3.15) 

Here all derivatives of (j)(t) are assumed to be expressed through 0(i) using (|3.13p . 

Let us simplify (|3.15p when = 4>max- At this point = 0, because 0(t) obeys an equation = —U'{cj)). Since the 
gap between energy and potential maximum is small and = \J E — U{cj)), so is also small. 

Substituting (jjmax = inD/sY/'^ into (|2.2ip . we obtain the value of Yi at (l)rnax'- Yi — 2nD + s{nD / s)^/'^ . Again, 
the first term is generated by the shear flow and at n ^ s/D dominates over the second, which is determined by 
fluctuations Using (j2.2ip and the fact that 4'{4)m) = we conclude that Yi = 0. 

Now we collect it all together: 

(5/2(0™) = 0, 

5/i(0m) = -T{2nD + s{nD/sY'^f. 

Thus, we've calculated the correction to the potential maximum. It ultimately determines correction to growth 
rate of (i?"). 

5\n = SE/D = 2nT{2nD + s{nD/sf/^f w ~2nT{nDs'^f'^ (3.16) 
Let us represent the answer in the relative form 

<5A„/A„ = -UinDs^f'^ = -^tX^. (3.17) 
3 9n 

Let us notice that correction p.l6p depends on n stronger than the main answer (|3.12p . However, ioi n s/D it 
remains much lesser the main answer, because we assume st <^ 1. 



IV. CONCLUSION 



We have studied the dynamics of inertial particles in a two-dimensional random flow with a strong shear component. 
We extracted information about the system from the statistical characteristics of the separation vector between two 
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particles R(i) . The moments of the distance grow exponentiaUy in time and the growth rates depend non-hnearly on 
n. We have calculated the first correction to the Lyapunov exponent, which is negative and of order tA^, where A is 
the inertialess Lyapunov exponent. This means that inertia reduces mixing and the relative correction to A is of order 
At. Let us note that for the shear-less situation, where inertialess Lyapunov exponent is 2D and the first correction 
due to inertia is Dt, relative correction is also of order Ar Q. 

We have also calculated A„ - the growth rates of the high moments of the inter-particle distance using a saddle-point 
approximation in the functional space (instanton method). The correction to the growth rate of the n-th moment 
due to inertia is found to be of order — nr(n£)s^)^/'^. Let us notice that the relative correction is proportional to 
SXn/Xn ^ W^rXn, SO the role of inertia reduces with n. 

Let us mention a few issues that are interesting to investigate further. 

As it is was shown in Q, in the isotropic model, as particles grow heavier, a phase transition to non-clustering 
phase may occur. We don't know yet if this transition exists in the presence of a strong shear. 

The velocity fluctuations we considered were white-correlated in time. Obviously this is an approximation, because 
any real noise is colored. It's unclear if a finite correlation time of the noise can lead to new qualitative effects. 

The Lyapunov exponents can be used to determine fractal characteristics of a cluster that is formed by particles, 
namely the so-called Lyapunov dimension (isl . [loj . For this purpose one has to calculate corrections to both Lyapunov 
exponents (in two dimensions) up to the second order in r, because the first order cancels in fractal dimension. 

In this paper we assumed that the parameters entering the model satisfy Dt <^ st <^ 1. Another limiting case, 
interesting for applications in atmospheric turbulence (though more complicated technically), is the case of larger 
inertia Dt <C 1 ^ sr. 
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